1  Análisis de expresión génica diferencial (RNA-Seq)

El análisis de la expresión génica diferencial (RNA-Seq) es una de las técnicas bioinformáticas que más auge ha experimentado en la última década, ya que permite cuantificar la expresión de un gen a nivel celular y comparar estas expresiones entre distintos grupos experimentales para identificar asociaciones entre genes y factores experimentales (como enfermedades o tratamientos).

El proceso de un análisis genético diferencial tiene distintas etapas, que van desde la obtención del las células a estudiar, pasando por la secuenciación del RNA que contienen, hasta el análisis estadístico de las frecuencias de expresión de cada gen.

Etapas del análisis de expresión génica diferencial

En este tutorial nos centraremos en el análisis estadístico de las frecuencias de expresión de los genes. Para ello partiremos de una matriz con las frecuencias o conteos de expresión de cada gen en el conjunto de muestras analizadas. Las filas de esta matriz corresponden a los genes observados y las columnas a las muestras analizadas.

	sample1	sample2	sample3	sample4	sample5	sample6
gene1	85	76	103	107	140	124
gene2	1	6	11	6	1	4
gene3	80	98	39	82	97	113
gene4	92	83	59	85	100	98
gene5	36	24	18	50	22	15
gene6	0	0	1	4	2	3

La mayor parte de los paquetes de Bioconductor que suelen encapsular esta matriz de frecuencias en un objeto del tipo SingleCellExperiment, que además de la matriz de frecuencias de expresión de los genes incluye también los grupos experimentales a los que pertenecen las muestras estudiadas y otra información que se va generando a medida que progresa el análisis.

El objeto SingleCellExperiment

El paquete SingleCellExperiment de Bioconductor define esta estructura de datos, y puede instalarse mediante el siguiente código

BiocManager::install('SingleCellExperiment')

Existen multitud de paquetes en el repositorio Bioconductor para realizar el análisis de expresión génica diferencial. En este tutorial veremos dos de los más usados: EdgeR y DesSeq2.

1.1 Análisis de expresión génica diferencial con EdgeR

1.1.1 Estructura de datos

El paquete EdgeR utiliza la estructura de datos DGEList para guardar la matriz de frecuencias de expresión génica. Para crear esta estructura se utiliza la función

• DGEList(frec, group = grupo): Crea una estructura del tipo DGEList con la matriz de frecuencias de expresión génica frec (con genes en filas y muestras en columnas) y el vector grupo con los grupos experimentales a los que pertenecen las muestras analizadas en las columnas de la matriz de frecuencias.

library(edgeR)
library(DEFormats)
frec <- simulateRnaSeqData()
grupo <- rep(c("A", "B"), each = 3)
dge <- DGEList(frec, group = grupo)
dge

An object of class "DGEList"
$counts
      sample1 sample2 sample3 sample4 sample5 sample6
gene1      85      76     103     107     140     124
gene2       1       6      11       6       1       4
gene3      80      98      39      82      97     113
gene4      92      83      59      85     100      98
gene5      36      24      18      50      22      15
995 more rows ...

$samples
        group lib.size norm.factors
sample1     A    42832            1
sample2     A    40511            1
sample3     A    39299            1
sample4     B    43451            1
sample5     B    40613            1
sample6     B    43662            1

Esta estructura de datos es una lista con dos atributos. El atributo counts contiene la matriz de frecuencias de expresión génica, y el atributo samples es un data frame información sobre las muestras estudiadas. Cada fila de este data frame se corresponde con una columna de la matriz de frecuencias y contiene las siguientes columnas

Columna	Descripción
group	Grupo experimental al que pertenece la muestra.
lib.size	tamaño de la librería (suma de frecuencias de la muestra).
norm.factors	Factor de normalización.

La estructura de datos DGEList puede contener opcionalmente el atributo genes con anotaciones de los genes observados en las filas de la matriz de frecuencias.

Para convertir esta estructura de datos en una del tipo DESeqDataSet se puede utilizar la función as.DESeqDataSet del paquete DEFormats.

dds <- as.DESeqDataSet(dge)
dds

class: DESeqDataSet 
dim: 1000 6 
metadata(1): version
assays(1): counts
rownames(1000): gene1 gene2 ... gene999 gene1000
rowData names(0):
colnames(6): sample1 sample2 ... sample5 sample6
colData names(3): group lib.size norm.factors

Ejemplo 1.1 Veamos cómo crear la estructura de datos correspondiente al estudio de Sheridan JM, et al. (2015) donde se obtuvieron datos de tres poblaciones de células (basal, progenitor luminal (LP) y luminal maduro (ML)) seleccionadas de las glándulas mamarias de ratones vírgenes hembra, cada una por triplicado. Los datos pueden obtenerse del repositorio Gene Expression Omnibus (GEO) mediante el identificador GSE63310.

En primer lugar descargamos los ficheros con las librerías de frecuencias.

# Descarga de datos
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="datos/GSE63310_RAW.tar", mode="wb") 
utils::untar("datos/GSE63310_RAW.tar", exdir = "datos")
files <- paste0("datos/", c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
  "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
  "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt"))
#for(i in paste0("datos/", files, ".gz"))
#  R.utils::gunzip(i, overwrite=TRUE)

A continuación creamos la estructura de datos DGEList.

# Creación de la estructura de datos DGEList
dge <- readDGE(files, columns=c(1,3))
dge

An object of class "DGEList"
$samples
                                                    files group lib.size
datos/GSM1545535_10_6_5_11 datos/GSM1545535_10_6_5_11.txt     1 32863052
datos/GSM1545536_9_6_5_11   datos/GSM1545536_9_6_5_11.txt     1 35335491
datos/GSM1545538_purep53     datos/GSM1545538_purep53.txt     1 57160817
datos/GSM1545539_JMS8-2       datos/GSM1545539_JMS8-2.txt     1 51368625
datos/GSM1545540_JMS8-3       datos/GSM1545540_JMS8-3.txt     1 75795034
datos/GSM1545541_JMS8-4       datos/GSM1545541_JMS8-4.txt     1 60517657
datos/GSM1545542_JMS8-5       datos/GSM1545542_JMS8-5.txt     1 55086324
datos/GSM1545544_JMS9-P7c   datos/GSM1545544_JMS9-P7c.txt     1 21311068
datos/GSM1545545_JMS9-P8c   datos/GSM1545545_JMS9-P8c.txt     1 19958838
                           norm.factors
datos/GSM1545535_10_6_5_11            1
datos/GSM1545536_9_6_5_11             1
datos/GSM1545538_purep53              1
datos/GSM1545539_JMS8-2               1
datos/GSM1545540_JMS8-3               1
datos/GSM1545541_JMS8-4               1
datos/GSM1545542_JMS8-5               1
datos/GSM1545544_JMS9-P7c             1
datos/GSM1545545_JMS9-P8c             1

$counts
           Samples
Tags        datos/GSM1545535_10_6_5_11 datos/GSM1545536_9_6_5_11
  497097                             1                         2
  100503874                          0                         0
  100038431                          0                         0
  19888                              0                         1
  20671                              1                         1
           Samples
Tags        datos/GSM1545538_purep53 datos/GSM1545539_JMS8-2
  497097                         342                     526
  100503874                        5                       6
  100038431                        0                       0
  19888                            0                       0
  20671                           76                      40
           Samples
Tags        datos/GSM1545540_JMS8-3 datos/GSM1545541_JMS8-4
  497097                          3                       3
  100503874                       0                       0
  100038431                       0                       0
  19888                          17                       2
  20671                          33                      14
           Samples
Tags        datos/GSM1545542_JMS8-5 datos/GSM1545544_JMS9-P7c
  497097                        535                         2
  100503874                       5                         0
  100038431                       1                         0
  19888                           0                         1
  20671                          98                        18
           Samples
Tags        datos/GSM1545545_JMS9-P8c
  497097                            0
  100503874                         0
  100038431                         0
  19888                             0
  20671                             8
27174 more rows ...

Añadimos la información del diseño experimental, en este caso el grupo experimental y el lote, al data frame de las muestras.

colnames(dge) <- substring(colnames(dge), 18, nchar(colnames(dge)))
# Añadimos el grupo experimental
dge$samples$group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP", 
                     "Basal", "ML", "LP"))
# Añadimos el lote
dge$samples$lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
dge

An object of class "DGEList"
$samples
                                   files group lib.size norm.factors lane
10_6_5_11 datos/GSM1545535_10_6_5_11.txt    LP 32863052            1 L004
9_6_5_11   datos/GSM1545536_9_6_5_11.txt    ML 35335491            1 L004
purep53     datos/GSM1545538_purep53.txt Basal 57160817            1 L004
JMS8-2       datos/GSM1545539_JMS8-2.txt Basal 51368625            1 L006
JMS8-3       datos/GSM1545540_JMS8-3.txt    ML 75795034            1 L006
JMS8-4       datos/GSM1545541_JMS8-4.txt    LP 60517657            1 L006
JMS8-5       datos/GSM1545542_JMS8-5.txt Basal 55086324            1 L006
JMS9-P7c   datos/GSM1545544_JMS9-P7c.txt    ML 21311068            1 L008
JMS9-P8c   datos/GSM1545545_JMS9-P8c.txt    LP 19958838            1 L008

$counts
           Samples
Tags        10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
  497097            1        2     342    526      3      3    535        2
  100503874         0        0       5      6      0      0      5        0
  100038431         0        0       0      0      0      0      1        0
  19888             0        1       0      0     17      2      0        1
  20671             1        1      76     40     33     14     98       18
           Samples
Tags        JMS9-P8c
  497097           0
  100503874        0
  100038431        0
  19888            0
  20671            8
27174 more rows ...

A continuación anotamos los genes de las filas de la matriz de frecuencias. Para ello debe utilizarse un paquete específico para el genoma de la especie de donde proviene el RNA (Mus.muculus en este caso para el genoma del ratón).

library(Mus.musculus)
geneid <- rownames(dge)
genes <- select(Mus.musculus, keys = geneid, columns = c("SYMBOL", "TXCHROM"), keytype = "ENTREZID")
# Eliminamos duplicidad de algunos genes
genes <- genes[!duplicated(genes$ENTREZID),]
dge$genes <- genes
dge

An object of class "DGEList"
$samples
                                   files group lib.size norm.factors lane
10_6_5_11 datos/GSM1545535_10_6_5_11.txt    LP 32863052            1 L004
9_6_5_11   datos/GSM1545536_9_6_5_11.txt    ML 35335491            1 L004
purep53     datos/GSM1545538_purep53.txt Basal 57160817            1 L004
JMS8-2       datos/GSM1545539_JMS8-2.txt Basal 51368625            1 L006
JMS8-3       datos/GSM1545540_JMS8-3.txt    ML 75795034            1 L006
JMS8-4       datos/GSM1545541_JMS8-4.txt    LP 60517657            1 L006
JMS8-5       datos/GSM1545542_JMS8-5.txt Basal 55086324            1 L006
JMS9-P7c   datos/GSM1545544_JMS9-P7c.txt    ML 21311068            1 L008
JMS9-P8c   datos/GSM1545545_JMS9-P8c.txt    LP 19958838            1 L008

$counts
           Samples
Tags        10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
  497097            1        2     342    526      3      3    535        2
  100503874         0        0       5      6      0      0      5        0
  100038431         0        0       0      0      0      0      1        0
  19888             0        1       0      0     17      2      0        1
  20671             1        1      76     40     33     14     98       18
           Samples
Tags        JMS9-P8c
  497097           0
  100503874        0
  100038431        0
  19888            0
  20671            8
27174 more rows ...

$genes
   ENTREZID  SYMBOL TXCHROM
1    497097    Xkr4    chr1
2 100503874 Gm19938    <NA>
3 100038431 Gm10568    <NA>
4     19888     Rp1    chr1
5     20671   Sox17    chr1
27174 more rows ...

1.1.2 Preprocesamiento

Una vez preparada la estructura de datos, el siguiente paso es el preprocesamiento de datos. Normalmente comprende las siguientes tareas:

1. Normalización de las frecuencias.
2. Eliminación de los genes con poca expresión.
3. Normalización de las distribuciones de frecuencias de expresión génicas.

1.1.2.1 Normalización de las frecuencias

El objetivo principal es normalizar las frecuencias de expresión génica para eliminar el efecto sobre las frecuencias de factores como la profundidad de secuenciado (a mayor profundidad de secuenciado mayores frecuencias) o el tamaño de los genes (a mayor tamaño mayores frecuencias). Generalmente se usan las siguientes transformaciones

• Frecuencias por millón (CPM). Se divide cada frecuencia por el tamaño en millones de la librería de la muestra a la que pertenece. Se realiza con la función cpm(dge).

• Logaritmo en base 2 de las frecuencias por millón (log2-CPM). Se calcula a partir de la anterior mediante la fórmula \(\log_2(CPM+\frac{2}{L})\), donde \(L\) es la longitud media de las librerías en millones. El término \(\frac{2}{L}\) que se añade permite calcular el logaritmo para frecuencias nulas (ya que el logaritmo de 0 no existe). Se realiza con la función cmp(dge, log = TRUE).

• Lecturas por kilobase de transcripción por millón (RPKM): Se realiza con la función rpkm(dge, longitud_genes) pasando la longitud de los genes.

• Fragmentos por kilobase de transcripción por millón (FPKM).

Las dos primeras no tienen en cuenta el tamaño de los genes, pero suelen usarse para la expresión génica diferencial donde el tamaño de los genes se supone constante en las distintas muestras.

Ejemplo 1.2 Siguiendo con el ejemplo anterior, calculamos las frecuencias por millón y el logaritmo de las frecuencias por millón.

cpm <- cpm(dge)
summary(cpm)

   10_6_5_11            9_6_5_11            purep53             JMS8-2        
 Min.   :    0.000   Min.   :    0.000   Min.   :   0.000   Min.   :   0.000  
 1st Qu.:    0.000   1st Qu.:    0.000   1st Qu.:   0.000   1st Qu.:   0.000  
 Median :    0.578   Median :    0.736   Median :   0.892   Median :   0.895  
 Mean   :   36.793   Mean   :   36.793   Mean   :  36.793   Mean   :  36.793  
 3rd Qu.:   19.536   3rd Qu.:   23.546   3rd Qu.:  24.221   3rd Qu.:  23.341  
 Max.   :27807.947   Max.   :11546.719   Max.   :7951.408   Max.   :7389.433  
     JMS8-3             JMS8-4              JMS8-5            JMS9-P7c       
 Min.   :   0.000   Min.   :    0.000   Min.   :   0.000   Min.   :   0.000  
 1st Qu.:   0.000   1st Qu.:    0.000   1st Qu.:   0.000   1st Qu.:   0.000  
 Median :   0.699   Median :    0.711   Median :   0.908   Median :   0.845  
 Mean   :  36.793   Mean   :   36.793   Mean   :  36.793   Mean   :  36.793  
 3rd Qu.:  23.827   3rd Qu.:   19.928   3rd Qu.:  21.439   3rd Qu.:  24.260  
 Max.   :7955.680   Max.   :29572.361   Max.   :9376.973   Max.   :7865.350  
    JMS9-P8c        
 Min.   :    0.000  
 1st Qu.:    0.000  
 Median :    0.752  
 Mean   :   36.793  
 3rd Qu.:   21.594  
 Max.   :16500.710  

lcpm <- cpm(dge, log = TRUE)
summary(lcpm)

   10_6_5_11          9_6_5_11          purep53             JMS8-2       
 Min.   :-4.5074   Min.   :-4.5074   Min.   :-4.50743   Min.   :-4.5074  
 1st Qu.:-4.5074   1st Qu.:-4.5074   1st Qu.:-4.50743   1st Qu.:-4.5074  
 Median :-0.6847   Median :-0.3589   Median :-0.09513   Median :-0.0901  
 Mean   : 0.1714   Mean   : 0.3312   Mean   : 0.43559   Mean   : 0.4089  
 3rd Qu.: 4.2913   3rd Qu.: 4.5601   3rd Qu.: 4.60081   3rd Qu.: 4.5475  
 Max.   :14.7632   Max.   :13.4952   Max.   :12.95700   Max.   :12.8513  
     JMS8-3            JMS8-4            JMS8-5            JMS9-P7c      
 Min.   :-4.5074   Min.   :-4.5074   Min.   :-4.50743   Min.   :-4.5074  
 1st Qu.:-4.5074   1st Qu.:-4.5074   1st Qu.:-4.50743   1st Qu.:-4.5074  
 Median :-0.4281   Median :-0.4064   Median :-0.07152   Median :-0.1704  
 Mean   : 0.3225   Mean   : 0.2529   Mean   : 0.40428   Mean   : 0.3708  
 3rd Qu.: 4.5772   3rd Qu.: 4.3199   3rd Qu.: 4.42513   3rd Qu.: 4.6031  
 Max.   :12.9578   Max.   :14.8520   Max.   :13.19491   Max.   :12.9413  
    JMS9-P8c      
 Min.   :-4.5074  
 1st Qu.:-4.5074  
 Median :-0.3300  
 Mean   : 0.2749  
 3rd Qu.: 4.4355  
 Max.   :14.0102  

1.1.2.2 Eliminación de genes con poca expresión

Aunque el objetivo del análisis de la expresión génica diferencial es detectar los genes que se expresan en un grupo experimental en comparación con los otros, cuando un gen no se expresa en ninguna de las muestras no tiene interés para el estudio y puede eliminarse.

Existen diferentes criterios para eliminar los genes con poca expresión. El paquete EdgeR incorpora la siguiente función

• filterByExpr(dge): Realiza un filtro de la estructura de datos dge los genes con poca expresión. Por defecto devuelve TRUE para cualquier gen con una frecuencia mayor o igual que 10 en al menos un número de muestras igual al del menor grupo experimental.

Ejemplo 1.3 Veamos cuántos genes no tienen expresión en ninguna muestra del ejemplo anterior.

# Número de genes con expresión nula.
table(rowSums(dge$counts) == 0)

FALSE  TRUE 
22026  5153 

El siguiente gráfico muestra la distribución del logaritmo en base 2 de las frecuencias por millón.

# Definimos una función para dibujar la distribución del logaritmo de las frecuencias por millón.
dist_logcpm <- function(lcpm) {
lcpm |> 
    as.tibble()  |>
    pivot_longer(everything(), names_to = "Muestra", values_to = "LogCPM") |>
    ggplot(aes(x = LogCPM, color = Muestra)) +
    geom_density() +
    labs(title = "Distribución del logaritmo en base 2 de las frecuencias por millón")
}

dist_logcpm(lcpm)

Warning: `as.tibble()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` instead.
ℹ The signature and semantics have changed, see `?as_tibble`.

Como se observa ha una porcentaje significativo de genes que tiene una expresión muy baja (valores negativos).

A continuación eliminamos de la estructura de datos los genes con poca expresión.

filtro <- filterByExpr(dge)
dge <- dge[filtro, , keep.lib.sizes = F]
dge

An object of class "DGEList"
$samples
                                   files group lib.size norm.factors lane
10_6_5_11 datos/GSM1545535_10_6_5_11.txt    LP 32857304            1 L004
9_6_5_11   datos/GSM1545536_9_6_5_11.txt    ML 35328624            1 L004
purep53     datos/GSM1545538_purep53.txt Basal 57147943            1 L004
JMS8-2       datos/GSM1545539_JMS8-2.txt Basal 51356800            1 L006
JMS8-3       datos/GSM1545540_JMS8-3.txt    ML 75782871            1 L006
JMS8-4       datos/GSM1545541_JMS8-4.txt    LP 60506774            1 L006
JMS8-5       datos/GSM1545542_JMS8-5.txt Basal 55073018            1 L006
JMS9-P7c   datos/GSM1545544_JMS9-P7c.txt    ML 21305254            1 L008
JMS9-P8c   datos/GSM1545545_JMS9-P8c.txt    LP 19955335            1 L008

$counts
        Samples
Tags     10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
  497097         1        2     342    526      3      3    535        2
  20671          1        1      76     40     33     14     98       18
  27395        431      771    1368   1268   1564    769    818      468
  18777        768     1722    2517   1923   3865   1888   1830     1246
  21399        810      977    2472   1870   2251   1716   1932      756
        Samples
Tags     JMS9-P8c
  497097        0
  20671         8
  27395       342
  18777       693
  21399       619
16619 more rows ...

$genes
  ENTREZID SYMBOL TXCHROM
1   497097   Xkr4    chr1
5    20671  Sox17    chr1
6    27395 Mrpl15    chr1
7    18777 Lypla1    chr1
9    21399  Tcea1    chr1
16619 more rows ...

lcpm <- cpm(dge, log = TRUE)
dist_logcpm(lcpm)

1.1.2.3 Normalización de las distribuciones de frecuencias de expresión génicas

Durante el proceso de secuenciación puede haber factores externos no biológicos que afecten a la expresión de los genes. Por ejemplo, la muestras del primer lote pueden tener mayores frecuencias que las del segundo lote. Como se supone que las muestras replicadas deben tener un rango y una distribución de frecuencias similares, en esta etapa se aplica otro procedimiento de normalización para garantizar que la distribución de frecuencias de cada muestra es similar. Para ello el paquete edgeR utiliza el método de las medias recortadas de los valores M por medio de la función

• calcNormFactors(dge, method = "TMM"): Calcular los factores de escalado del tamaño de las librerías de cada muestra. Estos factores se guardan automáticamente en el data frame con la información de las muestras dge$samples.

Ejemplo 1.4 Siguiendo con el mismo ejemplo, calculamos los factores de escalado para cada muestra.

dge <- calcNormFactors(dge, method = "TMM")
dge$samples$norm.factors

[1] 0.8943956 1.0250186 1.0459005 1.0458455 1.0162707 0.9217132 0.9961959
[8] 1.0861026 0.9839203

A continuación se muestran los diagramas de cajas de las distribuciones normalizadas tras aplicar los factores de escalado.

box_logcpm <- function(lcpm){
lcpm |> 
    as.tibble()  |>
    pivot_longer(everything(), names_to = "Muestra", values_to = "LogCPM") |>
    ggplot(aes(x = Muestra, y = LogCPM, fill = Muestra)) +
    geom_boxplot() +
    labs(title = "Distribución del logaritmo en base 2 de las frecuencias por millón")
}

box_logcpm(lcpm)

Como se puede apreciar, después de la normalización, todas las muestras presentan una distribución de frecuencias similar.

1.1.3 Análisis exploratorio

Una vez preprocesados los datos comienza el análiss estadístico propiamente dicho. En una primera fase se suele realizar un análisis exploratorio de los datos que suele incluir los siguientes análisis:

• Escalado multidimensional (análisis de componentes principales).

1.1.3.1 Escalado multidimensional

El escalado multidimensional mediante componente principales permite ver qué muestras son similares en cuando a la distribución de frecuencias de expresión génica. Los componentes principales son combinaciones lineales de los genes de la matriz de frecuencias que son ortogonales entre sí. El primer componente principal recoge la dimensión con mayor variabilidad de las frecuencias. El segundo recoge la dimensión don la mayor variabilidad no explicada por el primer componente principal y así sucesivamente. Normalmente los dos primeros componentes principales suelen recoger un porcentaje bastante alto de la variabilidad de las frecuencias. Al representar las muestras en estas dos dimensiones, las muestras más próximas entre sí, presentan una distribución de frecuencias de expresión génica similar. Para realizar esta representación se puede utilizar la siguiente función del paquete limma:

• plotMDS(dge): Realiza un diagrama de componentes principales de la matriz de frecuencias de la estructura de datosdge`.

Alternativamente, se pueden calcular los componentes principales mediante la función prcomp del paquete stats y luego usar la función autoplot del paquete ggfortify para dibujar el diagrama de componentes principales.

Ejemplo 1.5 A continuación se muestra cómo obtener el diagrama de componentes principales de la matriz de frecuencias de nuestro ejemplo.

plotMDS(dge, col = as.numeric(dge$samples$group), main = "Diagrama de componentes principales del logaritmo en base 2 de las frecuencias por millón")

library(ggfortify)
autoplot(prcomp(t(lcpm)), data = dge$samples, color = "group", loadings = TRUE, loadings.label = TRUE) +
labs(title = "Diagrama de componentes principales del logaritmo en base 2 de las frecuencias por millón")

Como se puede apreciar las muestras de cada grupo experimental aparecen agrupadas. La mayor diferencia (a lo largo del primer componente principal) se da entre el grupo basal y los otros dos grupos. Esto indica que cuando se realice el análisis de expresión génica diferencial habrá bastantes genes con diferencias de expresión significativa entre el grupo basal y los otros dos, mientras que no habrá tantos al comparar los grupos LP y ML.

Otra opción muy interesante es el paquete Glimma que permite dibujar un diagrama de componentes principales interactivo mediante la función

• glMDSPlot(lcpm): Dibuja un diagrama de componentes principales interactivo de la matriz de frecuencias lcpm.

library(Glimma)
#glMDSPlot(lcpm, groups = dge$samples[,c(2,5)])
dds <- DESeqDataSetFromMatrix(
  countData = dge$counts,
  colData = dge$samples,
  rowData = dge$genes,
  design = ~group
)

converting counts to integer mode

glimmaMDS(dds)

1.1.4 Análisis de expresión génica diferencial

La última etapa del análisis es la detección de los genes con una expresión significativamente diferente entre los grupos experimentales. Para ello se suelen utilizar modelos lineales o modelos lineales generalizados. En general, la estimación de los parámetros del modelo ajustado depende de la distribución teórica usada para modelizar las frecuencias de expresión génica. El paquete limma, por ejemplo, realiza un ajuste de modelo lineal suponiendo que las distribuciones de las variables implicadas son normales, mientras que el paquete edgeR modeliza las frecuencias de expresión de los genes observadas mediante la distribución binomial negativa.

El primer paso es definir la matriz del diseño del experimento, que define los grupos experimentales a comparar.

diseño <- model.matrix(~ 0 + dge$samples$group + dge$samples$lane)
colnames(diseño) <- gsub("dge\\$samples\\$group", "", colnames(diseño))
colnames(diseño) <- gsub("dge\\$samples\\$", "", colnames(diseño))
diseño

  Basal LP ML laneL006 laneL008
1     0  1  0        0        0
2     0  0  1        0        0
3     1  0  0        0        0
4     1  0  0        1        0
5     0  0  1        1        0
6     0  1  0        1        0
7     1  0  0        1        0
8     0  0  1        0        1
9     0  1  0        0        1
attr(,"assign")
[1] 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$`dge$samples$group`
[1] "contr.treatment"

attr(,"contrasts")$`dge$samples$lane`
[1] "contr.treatment"

Para definir los contrastes de ajuste por pares, el paquete limma utiliza la función

• makeContrast(pares)

contr.matrix <- makeContrasts(
   BasalvsLP = Basal-LP, 
   BasalvsML = Basal - ML, 
   LPvsML = LP - ML, 
   levels = colnames(diseño))
contr.matrix

          Contrasts
Levels     BasalvsLP BasalvsML LPvsML
  Basal            1         1      0
  LP              -1         0      1
  ML               0        -1     -1
  laneL006         0         0      0
  laneL008         0         0      0

En general en las distribuciones de frecuencias de expresión génica, se ha observado que la varianza no es independiente de la media. Los métodos que modelizan las frecuencias mediante el modelo binomial negativo asumen una relación cuadrática entre la media y la varianza. Con el paquete limma el ajuste lineal se realiza sobre el logaritmo en base 2 de las frecuencias por millón (log-CPM) que se suponen tienen una distribución normal. La relación entre la media y la varianza se realiza automáticamente mediante la función voom.

par(mfrow=c(1,2))
v <- voom(dge, diseño, plot=TRUE)
v

An object of class "EList"
$genes
  ENTREZID SYMBOL TXCHROM
1   497097   Xkr4    chr1
5    20671  Sox17    chr1
6    27395 Mrpl15    chr1
7    18777 Lypla1    chr1
9    21399  Tcea1    chr1
16619 more rows ...

$targets
                                   files group lib.size norm.factors lane
10_6_5_11 datos/GSM1545535_10_6_5_11.txt    LP 29387429    0.8943956 L004
9_6_5_11   datos/GSM1545536_9_6_5_11.txt    ML 36212498    1.0250186 L004
purep53     datos/GSM1545538_purep53.txt Basal 59771061    1.0459005 L004
JMS8-2       datos/GSM1545539_JMS8-2.txt Basal 53711278    1.0458455 L006
JMS8-3       datos/GSM1545540_JMS8-3.txt    ML 77015912    1.0162707 L006
JMS8-4       datos/GSM1545541_JMS8-4.txt    LP 55769890    0.9217132 L006
JMS8-5       datos/GSM1545542_JMS8-5.txt Basal 54863512    0.9961959 L006
JMS9-P7c   datos/GSM1545544_JMS9-P7c.txt    ML 23139691    1.0861026 L008
JMS9-P8c   datos/GSM1545545_JMS9-P8c.txt    LP 19634459    0.9839203 L008

$E
        Samples
Tags     10_6_5_11  9_6_5_11   purep53     JMS8-2    JMS8-3    JMS8-4    JMS8-5
  497097 -4.292165 -3.856488 2.5185849  3.2931366 -4.459730 -3.994060 3.2869677
  20671  -4.292165 -4.593453 0.3560126 -0.4073032 -1.200995 -1.943434 0.8442767
  27395   3.876089  4.413107 4.5170045  4.5617546  4.344401  3.786363 3.8990635
  18777   4.708774  5.571872 5.3964008  5.1623650  5.649355  5.081611 5.0602470
  21399   4.785541  4.754537 5.3703795  5.1220551  4.869586  4.943840 5.1384776
        Samples
Tags       JMS9-P7c  JMS9-P8c
  497097 -3.2103696 -5.295316
  20671  -0.3228444 -1.207853
  27395   4.3396075  4.124644
  18777   5.7513694  5.142436
  21399   5.0308985  4.979644
16619 more rows ...

$weights
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
[1,]  1.079413  1.332986 19.826915 20.273253  1.993686  1.395853 20.494977
[2,]  1.170357  1.456380  4.804866  8.660025  3.612508  2.626870  8.760149
[3,] 20.219073 25.573792 30.434759 28.528310 31.352260 25.743247 28.722497
[4,] 26.947557 32.505933 33.583128 33.232125 34.231754 32.354158 33.334340
[5,] 26.610864 28.501638 33.645479 33.206374 33.573492 31.996623 33.308490
          [,8]      [,9]
[1,]  1.107780  1.079413
[2,]  3.211473  2.541942
[3,] 21.200072 16.657930
[4,] 30.348630 24.259801
[5,] 25.171513 23.573305
16619 more rows ...

$design
  Basal LP ML laneL006 laneL008
1     0  1  0        0        0
2     0  0  1        0        0
3     1  0  0        0        0
4     1  0  0        1        0
5     0  0  1        1        0
6     0  1  0        1        0
7     1  0  0        1        0
8     0  0  1        0        1
9     0  1  0        0        1
attr(,"assign")
[1] 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$`dge$samples$group`
[1] "contr.treatment"

attr(,"contrasts")$`dge$samples$lane`
[1] "contr.treatment"

vfit <- lmFit(v, diseño)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")

A partir del modelo estimado se pueden obtener los genes subexpresados y sobreexpresados mediante el p-valor ajustado, tomando un nivel de significación del 5%.

summary(decideTests(efit))

       BasalvsLP BasalvsML LPvsML
Down        4648      4927   3135
NotSig      7113      7026  10972
Up          4863      4671   2517

Para la comparación entre los niveles de expresión en basal y LP, se encontró que 4648 genes están regulados negativamente en basal en relación con LP y 4.863 genes están regulados al alza en basal en relación con LP, un total de 9.511 genes significativamente diferenciados. Por otro lado, hay un total de 9598 genes significativamente diferenciados entre basal y ML (4927 genes regulados negativamente y 4671 regulados positivamente), y un total de 5652 genes significativamente diferenciados entre LP y ML (3135 regulados negativamente y 2517 regulados positivamente).

En algunos análisis se puede ser más conservador utilizando el log-fold-change (log-FC).

tfit <- treat(vfit, lfc=1)
dt <- decideTests(tfit)
summary(dt)

       BasalvsLP BasalvsML LPvsML
Down        1632      1777    224
NotSig     12976     12790  16210
Up          2016      2057    190

Para terminar se pueden mostrar los genes con diferencias más significativas.

basal.vs.lp <- topTreat(tfit, coef=1, n=Inf)
head(basal.vs.lp)

       ENTREZID SYMBOL TXCHROM     logFC  AveExpr         t      P.Value
12759     12759    Clu   chr14 -5.455444 8.856581 -33.55508 1.723731e-10
53624     53624  Cldn7   chr11 -5.527356 6.295437 -31.97515 2.576972e-10
242505   242505  Rasef    chr4 -5.935249 5.118259 -31.33407 3.081544e-10
67451     67451   Pkp2   chr16 -5.738665 4.419170 -29.85616 4.575686e-10
228543   228543   Rhov    chr2 -6.264208 5.485314 -29.07484 5.782844e-10
70350     70350  Basp1   chr15 -6.084738 5.247023 -28.26649 7.267694e-10
          adj.P.Val
12759  1.707586e-06
53624  1.707586e-06
242505 1.707586e-06
67451  1.739242e-06
228543 1.739242e-06
70350  1.739242e-06

basal.vs.ml <- topTreat(tfit, coef=2, n=Inf)
head(basal.vs.ml)

       ENTREZID  SYMBOL TXCHROM     logFC  AveExpr         t      P.Value
242505   242505   Rasef    chr4 -6.534668 5.118259 -35.08270 1.226399e-10
53624     53624   Cldn7   chr11 -5.495888 6.295437 -31.68918 2.774140e-10
12521     12521    Cd82    chr2 -4.690834 7.069637 -31.43063 2.913572e-10
20661     20661   Sort1    chr3 -4.931660 6.704459 -30.70113 3.558720e-10
71740     71740 Nectin4    chr1 -5.581283 5.164967 -30.59775 3.718157e-10
12759     12759     Clu   chr14 -4.686826 8.856581 -27.95760 7.687544e-10
          adj.P.Val
242505 1.236213e-06
53624  1.236213e-06
12521  1.236213e-06
20661  1.236213e-06
71740  1.236213e-06
12759  1.480434e-06

glimmaMA(tfit, coef=1, status=dt, main=colnames(tfit)[1],
         side.main="ENTREZID", counts=lcpm, groups=dge$samples$group)

External counts supplied using counts argument will be transformed to log-cpm by default. Specify transform.counts='none' to override transformation.

Warning in buildXYData(table, status, main, display.columns, anno, counts, :
count transform requested but not all count values are integers.

El mismo ajuste se puede realizar con el paquete edgeR.

dge <- estimateDisp(dge, diseño)
gfit <- glmFit(dge, diseño)
glrt <- glmLRT(gfit, diseño, contrast = contr.matrix)

#glimmaMA(glrt, dge = dge)

Y el volcano plot

glimmaVolcano(efit, dge = dge)

1.2 Análisis de expresión génica diferencial con DesSeq2

El paquete DESeq2 utiliza la estructura de datos DESeqDataSet para guardar la matriz de frecuencias de expresión génica. Para crear esta estructura se pueden utilizar las siguientes funciones

• DESeqDataSetFromMatrix(countData = frec, colData = grupo): Crea una estructura del tipo DESeqDataSet con la matriz de frecuencias de expresión génica frec (con genes en filas y muestras en columnas), y el data frame grupo con los grupos experimentales a los que pertenecen las muestras analizadas en las columnas de la matriz de frecuencias.

#library(DESeq2)
#dds <- DESeqDataSetFromMatrix(countData = frec, colData = dge$samples$group, design = ~ group)
#dds

1.3 Referencias

• Orchestrating Single-Cell Analysis with Bioconductor. (s. f.). Recuperado 12 de junio de 2023, de https://github.com/LTLA/OSCA

• Law, Charity, et al. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR.